Mon. Not. R. Astron. Soc. 0Q0.mi24l(2Q09) Printed 3 February 2010 (MN LffeX style file v2.2) 



The Orbital Evolution Induced by Baryonic Condensation in 
Triaxial Halos 



o 

(N 



Monica Valluri^*, Victor P. Debattista^, Thomas Quinn^, Ben Moore^ 

^ Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA 

^ RCUK Fellow: Jeremiah Horrocks Institute, University of Central Lancashire, Preston, PRI 2HE, UK 

Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195-1580, USA 
* Department of Theoretical Physics, University of Zurich, Winterthurerstrasse 190, CH-8057, Zurich, Switzerland 



Accepted - Dec, 11, 2009 



o 
u 

6 

(N 
> 

00 

O 

o 



X 



1 INTRODUCTION 



ABSTRACT 

Using spectral methods, we analyse the orbital structure of prolate/triaxial dark matter 
(DM) halos in iV-body simulations in an effort to understand the physical processes that drive 
the evolution of shapes of dark matter halos and elliptical galaxies in which central masses 
are grown. A longstanding issue is whether the change in the shapes of DM halos is the re- 
sult of chaotic scattering of the major family of box orbits that serves as the back-bone of a 
triaxial system, or whether they change shape adiabatically in response to the evolving galac- 
tic potential. We use the characteristic orbital frequencies to classify orbits into major orbital 
families, to quantify orbital shapes, and to identify resonant orbits and chaotic orbits. The use 
of a frequency-based method for distinguishing between regular and chaotic A^-body orbits 
overcomes the limitations of Lyapunov exponents which are sensitive to numerical discrete- 
ness effects. We show that regardless of the distribution of the baryonic component, the shape 
of a DM halo changes primarily due to changes in the shapes of individual orbits within a 
given family. Orbits with small pericentric radii are more likely to change both their orbital 
type and shape than orbits with large pericentric radii. Whether the evolution is regular (and 
reversible) or chaotic (and irreversible), depends primarily on the radial distribution of the 
baryonic component. The growth of an extended baryonic component of any shape results in 
aregular and reversible change in orbital populations and shapes, features that are not expected 
for chaotic evolution. In contrast the growth of a massive and compact central component re- 
sults in chaotic scattering of a significant fraction of both box and long-axis tube orbits, even 
those with pericenter distances much larger than the size of the central component. Frequency 
maps show that the growth of a disk causes a significant fraction of halo particles to become 
trapped by major global orbital resonances. We find that despite the fact that shape of a DM 
halo is always quite oblate following the growth of a central baryonic component, a significant 
fraction of its orbit population has the characteristics of its triaxial or prolate progenitor. 



The condensation of baryons to the centres of dark matter 
halo s is known to make them more spherical or axisymmet- 
ric toebattista et alj l2008l hereafter DOS). DOS found that the 
halo shape changes by A{b/a) > 0.2 out to at least half 
the virial radius. This shape change reconciles the strongly 
prolate-triaxial shapes found in coUisionless A'^-body simula- 



tions of the hierarchal growth of halos 
iBames & EfstathioJl987tlFreDk et alj 19881: 



Bardeenet al] 1 19861: 



Dubinski & Carlbera 



1991; Jing & Sutdl2002l : iBailin & Steimnetj 12005: Allgood etal] 
2006) with observation s, which g enerally find much round er halos 
(^Schw e i zer et al. 19S3i; Sackett & Sparke 1990; Franx & de~Zeeuwl 
19921: iHuizinga & van Albadal Il992l; iKuijken & Tremaind 1199? 



Franx etalj 1 1994 iBuote & Canizaresl [T994 iBartelmann et all 
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1995 


; iKochanek 19951: Ollind 1995L 19961; Schoenmakers et al. 


1997 


; Kooomans et al. 1998: lOlline & Merrifieldl 


200C: 


Andersen et al.' 2001'; 'Suote et al.' '2002'; Oauri et alj 


2003: 


Barnes & Sellwood 2003; Debattista 2003; lodice et al] 


2003; 


Diehl & Statlerll2007l;lBaneriee & Jodl2008l). 



What is the physical mechanism driving shape change? Op- 
tions suggested in the literature include two possibilities. The first is 
that the presence of a central mass concentration scatters box orbits 
that ser ve as the backbone of a t riaxial potential, render ing them 
chaotic dOerhard & Binnevlll985l ; lMerritt & Vallurilll996h . Chaotic 
orbits in a stationary potential do not conserve any integrals of 
motion other than the energy E and consequently are free to uni- 
formly fill their allowed equipotential surface. Since the potential 
is, in general, rounder than the density distribution chaotic diffu- 
sion results in evolution to a more oblate or ev en a spherical shape 
jMerritt & Ouinladl 19981 ; iKala potharakos 2008). The second pos- 
sibility is that the change of the central potential occurs because 
the growth of the baryonic component causes orbits of coUisionless 
particles in the halo to respond by chan ging their shapes in a reg- 
ular ( and therefore reversible) manner jHollev-Bockelmann et al.l 
l2002h . 

Time depende nce in a potential is al s o believed to result 
in ch aotic mixing dTerzic & Kandrupl |2004| ; iKandrup & Novotnvl 
I2OO4 !) and has been invoked as the mechanism that drives violent 
relaxation. However, a more recent analysis of mixing during a ma- 
jor merger showed that the rate and degree of mixing in energy 
and angular momentum are not consistent with chaotic mixing, but 
rather that particles retain strong memory of their initial energies 
and angular mome nta even in strongly time dependent potentials 
jValluri et al]|2007l) . 

One of the principal features of chaotic evolution is ir- 
reversibility. This irreversibility arises from two properties of 
chaotic orbits. First, chaotic orbits are exponentially sensitive to 
small changes in initial conditions even in a coUisionless sys- 
tem. Second, chaotic systems display the property of chaotic mix- 
ing (Lichtenberg & Lieberman 1992). Using this principle of ir- 
reversibility, DOS argued that if chaotic evolution is the primary 
driver of shape change, then if, subsequently, the central mass con- 
centration is artificially "evaporated", the system would not be able 
to revert to its original triaxial distribution. DOS showed that grow- 
ing baryonic components inside prolate/triaxial halos led to a large 
change in the shape of the halo. Despite these large changes, by ar- 
tificially evaporating the baryons, they showed that the underlying 
halo phase space distribution is not grossly altered unless the bary- 
onic component is too massive or centrally concentrated, or trans- 
fers significant angular momentum to the halo. This led them to 
argue that chaotic evolution alone cannot explain the shape change 
since such a process is irreversible. They speculated that at most 
only slowly diffusive chaos occurred in their simulations. Using test 
particle orbit integrations they also showed that box orbits largely 
become deformed, possibly changing into tube orbits, during disk 
growth, but do not become strongly chaotic. 



DOS employed irreversibility as a convenient proxy for the 
presence of chaos. In this paper we undertake an orbital analysis of 
some of the models studied by DOS to better understand the mech- 
anism that drives shape change. Our goal is to understand whether 
chaotic orbits are an important driver of shape change and if so 
under what conditions they are important. We also wish to under- 
stand how the orbital populations in halos change when a centrally 
concentrated baryonic component grows inside a triaxial dark mat- 
ter halo. Finally we would like to understand under what circum- 
stances orbits change their classification. 

This paper is organised as follows. In §|2]we describe the sim- 
ulations used in this paper and briefly describe three models from 
DOS as well as two additional simulations. In §[3] we describe the 
principal technique: Numerical Analysis of Fundamental Frequen- 
cies (NAFF) that we use to obtain frequency spectra and funda- 
mental frequencies and describe how these frequencies are used to 
characterise orbits. In §|4]we describe the results of our analysis 
of five different simulations. In §|5]we summarise our results and 
discuss their implications. 



2 NUMERICAL SIMULATIONS 

We formed prolate/triaxi al halos via mergers of systems, as 
described in iMoore et ^ ( l2004h . The initially spherical NFW 
jNavarro et alJl99 A halos were gen erated from a distribution func- 
tion using the method described in lKazantzidis et aL I ( l2004h with 
each halo composed of two mass species arranged on shells. The 
outer shell has more massiv e particles than the inner one, similar 
to the method described by IZemp et alj fcOOSi) . which allows for 
higher mass resolution at small radii. Our model halo A was gen- 
erated by the head-on merger of two prolate halos, themselves the 
product of a binary merger of spherical systems. The first merger 
placed the concentration c = 10 halos 800 kpc apart approach- 
ing each other at 50 kms^^, while the second merger starts with 
the remnant at rest, 400 kpc from an identical copy. The result- 
ing halo is highly prolate with a mild triaxiality. Halo model B 
was produced by the merger of two spherical halos starting at 
rest, SOO kpc apart and is prolate, with {b/a) = (c/a) ~ 0.58. 
Halo A has lb/a) ~ 0.45 and (c/a) ~ 0.35 while halo B has 
(b/a) = (c/a) ~ 0.58 (see Figure 3 of DOS for more details). Both 
halos A and B consist of 4 x 10*' particles. The outer particles are 
~ 18 times more massive in halo A and ~ 5 times more massive 
in halo B. A large part of the segregation by particle mass persists 
after the mergers a nd the small ra dius regions are dominated by low 
mass particles (cf. 'Dehnen''2005). We used a softening parameter 
e = 0.1 kpc for all halo particles. The radius, r2oo, at which the 
halo density is 200 times the mean density of the Universe and the 
total mass within this radius, M'zoo, are given in Table 1. 

Once we produced the prolate/triaxial halos, we inserted a 
baryonic component, either a disk of particles that remains rigid 
throughout the experiments or softened point particles. The param- 
eters that describe the distribution of the baryonic components are 
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Run Number 


Run 


Halo 


r200 


M200 




fb 


Rb 






(from DOS) 


Description 




[kpc] 




[1O"M0] 




[kpc] 


[Gyr] 


[Gyr] 


SAl 


Triax+Disk 


A 


215 


4.5 


1.75 


0.039 


3.0 


5 


2.5 


P,A3 


Triax+Bulg 


A 


215 


4.5 


1.75 


0.039 


1.0 


5 


2.5 


P/B2 


Proh+Ellip 


B 


106 


0.65 


0.7 


0.108 


3.0 


10 


4 


P;A4 


Triax+hardpt 


A 


215 


4.5 


1.75 


0.039 


0.1 


5 


2.5 


P;B3 


Prolt+hardpt 


B 


106 


0.65 


0.35 


0.054 


0.1 


5 


5 



Table 1. The simulations in this paper. Mt is the mass in baryons and ft is the baryonic mass fraction. For the particle simulations (PjB2, PiB3, P; A3, P;A4), 
refers to the softening of the spherical baryonic distribution particle(s). For simulation SAl, refers to the scale length of the baryonic disk. 



given in Tablel. In four of the models (PiA3, P;A4, P/B2 and 
PiB3) tlie baryonic component is simply a softened point mass witli 
softening scale length given by Rt- In model SAl, the density dis- 
tribution of the disis was exponential with scale length of the bary- 
onic component i?b and Gaussian scale-height Zb/^b ~ 0.05. The 
disk was placed with its symmetry axis along the triaxial halo's 
short axis in model SAl (additional orientations of the disk relative 
to the principal axes were also simulated but their discussion is de- 
ferred to a future paper). Initially, the disk has negligible mass, but 
it grows adiabatically and linearly with time to a mass Alt during 
a time tg. After this time, we slowly evaporated it during a time te- 
We stress that this evaporation is a numerical convenience for test- 
ing the effect of chaos on the system, and should not be mistaken 
for a physical evolution. The disk is composed of 300A' equal-mass 
particles each with a softening e = 100 pc. From t = lo tg + 
the halo particles are free to move and achieve equilibrium with the 
baryons as their mass changes, but all disk particles are frozen in 
place. The masses of models with single softened particles are also 
grown in the same way; these are models P/B2, PiB3 from DOS, 
and P;A3 and P;A4 which are new to this paper. DOS's naming 
convention for these experiments used "P" subscripted by "f" for 
particles frozen in place and by "1" for live particles free to move. 

Three different baryonic components are grown in the triax- 
ial halo A: in model SAl the baryons are in the form of a disk 
grown perpendicular to the short axis and the model is referred to as 
Triax+Disk; in model PiA3 the baryonic component is a softened 
central point mass resembling a bulge and the model is referred to 
as Triax+Bulg; finally in model P;A4 the baryonic component is a 
hard central point mass with a softening of 0. 1 kpc and the model 
is referred to as Triax+hardpt. Two different baryonic components 
are grown in the prolate halo B: in model P/B2 the baryonic com- 
ponent loosely resembles an elliptical galaxy so this run is referred 
to as Prolt+Ellip; in model P;B3 the baryonic component is a hard 
central point mass with a softening of 0. 1 kpc and is referred to as 
Prolt+hardpt. For model P/B2, DOS showed that there is no sig- 
nificant difference in the evolution if the central particle is live in- 
stead of frozen, all other things being equal. P; A3 was constructed 
specifically for this paper in order to have a triaxial halo model with 
a moderately soft spherical baryonic distribution which can be con- 
trasted with the prolate halo model P/B2, while P;A4 is a triaxial 



halo model which can be contrasted with the prolate halo model 
PiB3. 

For each model there were 5 phases in evolution. The initial 
triaxial or prolate halo without the baryonic component is referred 
to as phase a. There is then a phase (of duration tg) during which 
the halo's shape is evolving as the baryonic component is grown 
adiabatically. We do not study orbits in this phase since the po- 
tential is evolving with time. When the baryonic component has 
finished growing to full strength, and the halo has settled to a new 
equilibrium the model is referred to as being in phase b. The bary- 
onic component is "adiabatically evaporated" over a timescale of 
duration t^ listed in Table [T] Again we do not study orbits dur- 
ing this period when the potential is evolving with time. After the 
baryonic component has been adiabatically evaporated completely 
and the halo has returned to an equilibrium configuration the halo 
is referred to as being in phase c. We only study the orbits in the 
halo during the three phases when the halo is in equilibrium and is 
stationary (not evolving with time). 

The growth of the baryonic component induces several 
changes in the distributions of the DIVI halo particles: first is an 
increase in central density relative to the original NFW halo due to 
increase in the depth of the central potential (an effect commonly 
referred to as "baryonic compression"); second, the halos become 
more oblate especially within 0.3r„i,,- The details of the changes in 
the density and velocity distributions of DIVI particles differ slightly 
depending on the nature of the baryonic component (DOS). 

All the simulations in this paper, which are listed in Table [T] 
were evol ved with PKD GRAV an efficient, multi-stepping, parallel 
tree code jStadelll200ll) . We used cell opening angle for the tree 
code of — 0.7 throughouQ. Additional details of the simulations 
can be found in DOS. 



2.1 Computing Orbits 

In each of the halos studied we selected a subsample of between 
1000-6000 particles and followed their orbits in each of the three 



^ Opening angle 9 is used in tree codes to deter mine how long-range forces 
from particles acting at a point are accumulated iBames &Hui 19861) . 
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stationary phases of the evolution described in the previous sec- 
tion. The particles were randomly chosen in the halos at t = such 
that they were inside a fixed outer radius (either 100 or 200 kpc). 
Since the particles were selected at random from the distribution 
function, they have the same overall distribution as the entire dis- 
tribution function within the outer radius selected. We integrated 
the motion of each a test particle while holding all the other parti- 
cles fixed in place. We used a fixed timestep of 0. 1 Myr and inte- 
grated for 50 Gyr, storing the phase space coordinates of each test 
particle every 1 Myr. We used such long integration times to en- 
sure we are able to obtain accurate measurements of frequencies 
(as described in the next section). We carried out this operation for 
the same subset of particles at phases a, b and c. In model SAl 
(Triax+Disk) we integrated the orbits of 6000 particles which in 
phase a were within r = 200 kpc. In model PjAS (Triax+Bulg) 
and P;A4 (Triax+hardpt) we considered a subsample of 5000 par- 
ticles starting within r — 100 kpc. In models P/B2 and P;B3 we 
considered orbits of 1000 particles within r = 200 kpc. We inte- 
grated their orbits as above but we used a smaller timestep St = 10^ 
years in the case of P!A4 and PiB3, which had harder central point 
masses. The orbit code computes forces in a frozen potential using 
an integration scheme that uses forces calculated from the PKD- 
GRAV tree; we used the orbit integration parameters identical to 
those used for the evolution of the self-consistent models. 



3 FREQUENCY ANALYSIS 

In a 3-dimensional galactic potential that is close to integrable, all 
orbits are quasi-periodic. If an orbit is quasi-periodic (or regular), 
then any of its coordinates can be described explicitly as a series. 



(1) 



where the oj^'s are the oscillation frequencies and the Ai^'s are the 
corresponding amplitudes. In a three dimensional potential, each 
cok can be written as an integer linear combination of three fun- 
damental frequencies LUi,iU2,ijJ3 (one for each degree of freedom). 
If each component of the motion of a particle in the system (e.g. 
x{t)) is followed for several (~ 100) dynamical times, a Fourier 
transform of the trajectory yields a spectrum with discrete peaks. 
The locations of the peaks in the spectrum correspond to the fre- 
quencies LUk and their amplitudes Ak can be us ed to compute 
the linearly independent fundamental frequencies jBoozeJ 1 1 9821 : 
iKuo-Petravic e t al.'igSsVBinn ev & Trema ine"2008\ 

Binnev & Spergel (1982. Il984l) applied this method to 



galactic potentials and obtained the frequency spectra using a 
least s quares technique to measure the frequencies ojk- iLaskaj 
( Il990lll993h developed a significantly improved numerical tech- 
nique (Numerical Analysis of Fundamental Frequencies, hereafter 
NAFF) to decompose a complex time series of the phase space tra- 
jectory of an orbit of the form x(t) + ivxit), (where is the 



velocity along the x coordinate). IValluri & Merrit 3 i 19981) devel- 
oped their own implementation of this algorithm that uses integer 
programming to obtain the fundamental frequencies from the fre- 
quency spectrum. In this paper we use this latter implementation of 
the NAFF method. 

We re fer readers to the above papers and to Section 3.7 of 
iBinnev & T remaine (2008) for a detailed discussion of the main 
idea behind the recovery of fundamental frequencies. For complete- 
ness we provide a brief summary here. The NAFF algorithm for fre- 
quency analysis allows one to quickly and accurately compute the 
fundamental frequencies that characterise the quasi-periodic mo- 
tion of regular orbits. The entire phase space at a given energy can 
then be represented by a frequency map which is a plot of ratios 
of the fundamental frequencies of motion. A frequency map is one 
of the easiest ways to identify families of orbits that correspond to 
resonances between the three degrees of freedom. 

The structure of phase space in 3-dimensional galactic poten- 
tials is quite complex and we summarize some of its properties here 
to enable the reader to more fully appreciate the results of the analy- 
sis that follows. When an integrable potential is perturbed, its phase 
space structure is altered, resultin g in the appearance of resonances 
jLichtenberg & LiebermanI 19921) . Resonances are regions of phase 
space where the three fundamental frequencies are not linearly in- 
dependent of each other, but two or more of them are related to 
each other via integer linear relations. As the perturbation in the 
potential increases, the potential deviates further and further from 
integrability, and a larger and larger fraction of the phase space be- 
comes associated with resonances. 

In a three dimensional potential, orbits that satisfy one reso- 
nance condition such as luj^ + mujy -\- nuiz ~ are referred to as 
"thin orbit" resonances since t hey cover the surface of a two dimen- 
sional surface in phase space jMerritt & Vallurill 19991) . If two inde- 
pendent resonance conditions between the fundamental frequencies 
exist, then the orbit is a closed periodic orbit. Orbits that have fre- 
quencies close to the resonant orbit frequencies are said to be res- 
onantly trapped. Such orbits tend to have properties similar to that 
of the parent resonance, but get "thicker" as their frequencies move 
away from the resonance. At the boundary of the region of phase 
space occupied by a resonant family is a region called the "separa- 
trix". The separatrix is the boundary separating orbits with different 
orbital characteristics. In this case it is the region between orbits 
that have frequencies that are similar to the resonant orbits and or- 
bits that are not resonant. Chaotic orbits often occur in a "stochastic 
layer" close to resonances and at the intersections of resonances. In 
fact one of the primary factors leading to an incr ease in the frac - 
tion of chaotic orbits is the overlap of resonances jChirikov|[T979h . 
Chaotic orbits that are close to a resonant family are referr ed to 
as "resonantly trapped" or "sticky orbits" l lHabib et al.ll 1997 1) and 
are often only weakly chaotic. Orbits that are "sticky" behave like 
the resonant parent orbit for extremely long times and therefore do 
not diffuse freely over their energy surface or undergo significant 
chaotic mixing. 



© 2009 RAS, MNRAS Q0Q.mt24l 



The Orbital Evolution Induced by Baryonic Condensation in Triaxial Halos 5 



The frequency analysis method allows one to map the phase 
space structure of a distribution function and to easily identify the 
most important resonances by plotting ratios of pairs of frequen- 
cies (e.g. ujx/uJz vs. ujy/ujz) for many thousands of orbits in the 
potential. In such a frequency map, resonances appear as straight 
lines. Stable resonances appear as filled lines about which many 
points cluster, and unstable resonances appear as "blank" or de- 
populated lines. The strength of the resonances can be determined 
by the number of orbits that are associated with them. 



3.1 Overcoming microchaos in A^-body simulations 

In TV-body systems like those considered in this paper, the galactic 
mass distribution is realised as a discrete set of point masses. The 
discretization of the potential is known to result in exponential devi- 
ation of near by orbits, even in systems where all orbits are expected 
to be regular (Miller 1964; Goodman et al. 1993; Kandrup & Smith 
[199 1 ; Valluri & Merritt 2000; Hemsendorf & Merritt 2002). How- 
ever as the number of particles in a simulation is increased, and 
when point masses are softened, the majority of orbits begin to ap- 
pear regular despite the fact that t heir non-zero Lyapunov expo- 
nent implies that they are chaotic. iHemsendorf & Merritll ilOO'j) 
showed that this Lyapunov exponent saturates at a finite value 
beyond a few hundred particles and corresponds to an e-folding 
timescale of 1/20 of a system crossing time (for systems with 
A'^ ~ 10^ particles). Despite having large Lyapunov exponents (i.e. 
short e-folding times) these orbits behave and look much like regu- 
lar orbits tKandrup & Sideris .2001; . Kandrup & Siopis 20031) . This 
property of A'^-body orbits to have n on-zero Lyapunov expone nts 
has been referred to as "microchaos" /Kandrup & Sideris''200!j ^ or 
the "M iller Instability" (Hemsendorf & Merritt 2002; Valluri e VsH 
and suggests that Lyapunov exponents, while useful in con- 
tinuous potentials, are not a good measure of chaotic behavior re- 
sulting from the global potential when applied to A'^-body systems. 
This is a strong motivation for our use of a frequency based method 
which, as we demonstrate, is extremely effective at distinguishing 
between regular and chaotic orbits and is apparently largely unaf- 
fected by microchaos. 

We now discuss how frequency analysis can be used to distin- 
guish between regular and chaotic orbits. In realistic galactic po- 
tentials most chaotic orbits are expected to be weakly chaotic and 
lie close to regular orbits mimicking their behaviour for long times. 
The rate at which weakly chaotic orbits c hange their orb ital fre- 
quencies can be used as a measure of chaos. lLaskaJh993l) showed 
that the change in the fundamental frequencies over two consecu- 
tive time intervals can be used as a measure of the stochasticity of 
an orbit. This method has be en used to study the phase spac e struc- 
ture in galactic potent ials dPapaphilippou & LaskaJ 1 19961 1 19981 ; 
IValluri & Merrittlll998l) . Examples of frequency spectra for each 
component of motion , and their resolution into 3 fun damental fre- 
quencies are given bv lPapaphilippou & Laskaj d 19981 for different 
types of orbits. For each time series the spectrum is analysed and 



the three fundamental frequencies are obtained. In Cartesian coor- 
dinates the frequencies would be ijjx,ijjy,u}z- 

For each orbit we therefore divide the integration time of 50 
Gyr time into two consecutive segments and use NAFF to compute 
the fundamental frequencies ujx,iLJy, tOz (note that all frequencies 
in this paper are in units of Gyr^^, therefore units are not explicitly 
specified everywhere). We compute the three fundamental frequen- 
cies ijJx{ti),LOy{ti),ujz{ti) and L0x{t2) , ujy{t2) , 002(12) in each of 
the two intervals ti and t2 respectively. We compute the "frequency 
drift" for each frequency component as: 

log(A/.)=log|^^^^M_^|, (2) 
log(A/„)=log|^;^^M_^|, (3) 

UJy{tl) 

log A/, = log ' ^ ' . (4) 

We define the frequency drift parameter log(A/) (logarithm to 
base 10) to be the value associated with the largest of the three 
frequencies fx, fy, fz ■ The larger the value of the frequency drift 
parameter, the more chaotic the orbit. 

Identifying truly chaotic behavior however also requires that 
we properly account for numerical noise. In previous studies orbits 
were integrated with high numerical precision for at least 100 or- 
bital periods, resulting in highly accurate frequency determination. 
For instance. IValluri & MerrittI l[l998) found that orbital frequen- 
cies in a triaxial potential could be recovered with an accuracy of 
10"^'^ for regular orbits and 10"'* — 10"® for stochastic orbits using 
integration times of at least 50 orbital periods per orbit. 

In order to use frequency analysis to characterise orbits as reg- 
ular or chaotic in A'^-body systems, it is necessary to assess the nu- 
merical accuracy of orbital frequencies obtained by the NAFF code. 
To quantify the magnitude of frequency drift that arises purely from 
discretization effects (the microchaos discussed above) we select 
a system that is spherically symmetric and in dynamical equilib- 
rium. All orbits in a smooth spherically sym metric potentia l are 
rosettes confined to a single plane dBinnev & Tremainell2008h and 
are regular. Hence any drift in orbital frequencies can be attributed 
entirely to discretization errors (including minute deviations of the 
A'^-body potential from perfect sphericity). As a test of our applica- 
tion of the NAFF code to A'^-body potentials we analyse orbits in 
spherical NFW halos of two different concentrations (c = 10 and 
c — 20). The halos are represented by lO'' particles and have mass 
~ 2 X lO*'^ Mq . Particles in both cases come in two species with 
softening of 0. 1 kpc and 0.5 kpc. We carried out the frequency anal- 
ysis of 1000 randomly selected orbits which were integrated for 50 
Gyr in the frozen A'^-body realisations of each of the NFW halos. 

Figure [T] shows the distribution of values of log( A/) for both 
halos. In both cases the distribution has a mean value of log( A/) — 
—2.29, with standard deviations of 0.58 (for the c = 10 halo) and 
0.54 (for the c = 20 halo). Both distributions are significantly 
skewed toward small values of log(A/) (skewness = -0.85) and 
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Figure 1. Distributions of frequency drift parameter log(A/) for 1000 or- 
bits in two different spherical NFW halos. Despite the difference in con- 
centration c = 10 and c = 20 the distributions are almost identical having 
a mean of log(A/) = —2.29 and a standard deviation a ~ 0.56, with a 
significant skewness toward small values of log(A/). 



are more peaky than Gaussian (kurtosis = 1.95). Despite the fact 
that the two NFW halos have different concentrations, the distribu- 
tions of log(A/) are almost identical, indicating that our chaotic 
measure is largely independent of the central concentration. 

To define a threshold value of log(A/) at which orbits are 
classified as chaotic we note that 99.5% of the orbits have values 
of log(A/) < —1.0. Since all orbits in a stationary spherical halo 
are expected to be regular, we attribute all larger values of log (A/) 
to numerical noise arising from the discretization of the potential. 
Henceforth, we classify an orbit in our A'^-body simulations to be 
regular if it has log(A/) < —1.0. 

To accurately measure the frequency of an orbit it is necessary 
to sample a significant part of its phase-space structure (i.e. the sur- 
face of a 2-torus in a spherical potential or the surface of the 3-torus 
in a triaxial potential). Valluri & Merritt ( 1998) showed that the ac- 
curacy of the frequency analysis decreases significantly when or- 
bits were integrated for less than 20 oscillation periods. Inaccurate 
frequency determination could result in misclassifying orbits as 
chaotic (since inaccurate frequency measurement can also lead to 
larger frequency drifts). We test the dependence of log( A/) on the 
number of orbital periods Up by plotting the frequency drift param- 
eter against the largest orbital frequency (for orbits with Up > 20) 
in both NFW halos in Figure Qy. We use lj instead of Up since 



We use the fractional change in the largest of the three fundamental fre- 
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Figure 2. log(A/) versus ui (in units of Gyr"'^) for orbits with rip > 20 
in the two spherical NFW halos. Stars are for the halo with c = 10 and the 
open circles are for the halo with c = 20. The 1000 particles are binned in 
ui so that each bin contains the same number of particles. The vertical error 
bars represent the standard deviation in each bin. The straight lines are fits 
to the data, the slopes of both lines are consistent with zero. 



rip (X oj but is harder to compute accurately. Particles are binned in 
equal intervals in oj and the error bars represent the standard devi- 
ation in each bin. The straight-lines are best fits to the data-points. 
The slopes of the correlation for the c = 10 halo (solid line) and 
for the c = 20 halo (dot-dashed line) are both consistent with zero, 
indicating that log(A/) is largely independent of u (and hence of 
Tip). Henceforth we only use orbits which execute more than 20 
orbital periods in the 50 Gyr over which they are integrated. The 
excluded orbits lie predominantly at large radii and are not signif- 
icantly influenced by the changes in the inner halo that are investi- 
gated here. This rejection criterion affects about 25% of the orbits 
in the triaxial dark matter halos that we consider later. 

The effect of central concentration on the accuracy of fre- 
quency estimation is of particular concern during phase b, when 
the potential is deepened due to the growth of a baryonic compo- 
nent. In this phase, frequencies of those orbits which are strongly 
influenced by the deepened potential are increased. Consequently 
some orbits execute many more orbital periods during phase b than 
they do in phase a or phase c. However we have fixed the orbital 
sampling time period (not integration timestep) to 1 Myr in all 
phases. In principle coarse time sampling should not be a concern 
since the long integration time can still ensure a proper coverage 



quencies measured over two contiguous time intervals (frequency drift) as 
a measure of chaos ( Laskar 1990). For situations where a large fraction of 
orbits is resonant, it may be more appropriate to use the smallest of the three 
frequencies or the component with the largest amplitude. 
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of the phase-space torus. To ensure that the samphng frequency 
per orbital period does not significantly alter the frequency esti- 
mation we re-simulated one model (SAl) in phase b and stored 
the orbits 5 times more often (i.e. at time intervals of 0.2 Myr). 
We compared the frequencies of orbits computed for the low (ujl) 
and high {ujh) time-sampling runs. Figure [3] (fe//) shows that there 
is a strong correspondence between frequencies obtained with the 
two different samplings. We also found (Pi^^middle) that the fre- 
quency drift parameter log(A/) obtained from the two runs are 
highly correlated although there is some increase in scatter for or- 
bits with values of log(A/) > —2. Since the scattered points lie 
roughly uniformly above and below the 1 : 1 correlation line, there 
is no evidence that the higher sampling rate gives more accurate 
frequencies. The right panel shows that the overall distribution of 
Iog( A/) is identical for the two runs. We find that 95% of the par- 
ticles showed a frequency difference < 0.1% between the two dif- 
ferent sampling rates. From these tests we conclude that our choice 
of sampling rate in the phase b is unlikely to significantly affect 
the frequency measurements of the majority of orbits. We therefore 
adopt the lower orbit sampling frequency for all the analysis that 
follows. 



In addition to classifying orbits into these three broad categories, 
they showed that if one or more of the fundamental frequencies 
is an integer linear combination of the other frequencies, the or- 
bit can be shown to be resonant (either a periodic orbit or an open 
resonance). We followed the scheme outlined by CA98 to develop 
our own algorithm to classify orbits as boxes, long-axis tubes (ab- 
breviated as L-tubes) and short-axis-tubes (abbreviated as S-tubes) 
and to also identify orbits that are associated with low-order reso- 
nances. We do not describe the classification scheme here since it is 
essentially identical to that described by CA98, the main difference 
lies in that we use NAFF to obtain the fundamental frequencies of 
orbits in the A^-body model, wh ereas they used a method based on 
that of lBinnev & Spergell jl984h . We tested our automated classifi- 
cation by visually classifying 60 orbits that were randomly selected 
from the different models. We then ran our automated orbit classi- 
fier on this sample, and compared our visual classification with that 
resulting from the automated classifier. The two methods agreed for 
58/60 orbits (a 96% accuracy rate assuming that the visual classi- 
fication is perfectly accurate). Hereafter we assume that our auto- 
mated classification is accurate 96% of the time and therefore any 
orbit fractions quoted have an error of ±4%. 



3.2 Orbit classification 

ICarpintero & Aguil^ ( Il998l , hereafter CA98) showed that once a 
frequency spectrum of an orbit is decomposed into its fundamental 
frequencies, the relationships between the values of the frequencies 
{uix , ujy , uiz) can be used to classify the orbits in a triaxial potential 
into the major orbit families as boxes, long (x) axis tubes and short 
(2) axis tubes. (CA98 point out that it is difficult to distinguish be- 
tween the inner long-axis tubes and outer long-axis tubes from their 
frequencies alone. Therefore we do not attempt to distinguish be- 
tween these two families with our automatic classification scheme.) 



3.3 Quantifying orbital shapes 

In any self-consistent potential the distribution of shapes of the ma- 
jority of the orbits match the overall shape of the density profile. 
The elongation along the major axis is provided either by box orbits 
or by inner L-tubes. The ratios of the fundamental frequencies of 
orbits can be used to characterise their overall shape. Consider a tri- 
axial potential in which the semi-major axis (along the x-axis) has 
a length ax, the semi-intermediate axis has length ay, and the semi- 
minor axis has a length a^. The fact that ax > ay > a^ implies that 
the oscillation frequencies along these axes are \ijJx \ < \ijJy \ < \ijJz\ 
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for any (non resonant) orbit with tiie same over-all shape as the 
density distribution (we consider only the absolute values of the 
frequencies since their signs only signify the sense of oscillation). 
We can use this property to define an average "orbit shape param- 
eter" (Xs) for any orbit. For an orbit whose overall shape matches 
the shape of the potential, 

\uj-^ \ > \ujy\ > \uJx\ 
_^ \0Jy\ ^ 1^x1 

xs = f^-H>0- (5) 

The orbit shape parameter Xa is positive for orbits with elongation 
along the figure. The larger the value of Xs, the greater the degree 
of elongation along the major axis. Very close to the centre of the 
potential it is possible for orbits to have greater extent along the y 
axis than along the x axis, as is sometimes the case with outer L- 
tubes. For such orbits Xs is slightly negative. An orbit for which all 
frequencies are almost equal would enclose a volume that is almost 
spherical. For such an orbit, Xs ~ (which we refer to as "round"). 
Note that orbits which are close to axisymmetric about the short 
{z) axis (i.e. the S-tubes) also have Xs ~ because Ux ~ ojy 
regardless of the value of ujz- Our definition of shape parameter 
does not permit us to distinguish between truly round orbits for 
which uix ~ ^ and S-tubes, but both contribute to a more 
oblate axisymmetric potential. 



4 RESULTS 

For every model in Table [T] the three fundamental frequencies 
UxyUiy, LUz of each of the orbits in a selected subsample were com- 
puted separately in each of the three phases. For each orbit the 
largest of the three fundamental frequencies is assumed to repre- 
sent the dominant frequency of motion. The absolute value of this 
quantity is referred to as the largest fundamental frequency: we use 
Ua, uJbyOJc to refer to the largest fundamental frequencies of an or- 
bit in each of the three phases a, b, c. In addition to computing the 
fundamental frequencies over the entire 50 Gyr interval, we split 
the interval into two equal halves and computed the frequencies in 
each to compute the frequency drift parameter log(A/) defined in 
^ All orbits with log (A/) < —1.0 are identified as regular and 
the rest are identified as chaotic. For each orbit we also compute 
the total energy {E), the absolute value of the total specific angu- 
lar momentum (|jtot |), the number of orbital periods (rip), and the 
pericenter and apocenter distance from the centre of the potential 

(^pcri, ^apo). 

In this section we consider results of five simulations, SAl 
(Triax+Disk), P,A3 (Triax+Bulg), P,A4 (Triax+hardpt), P/B2 
(Frolt+Ellip), and P;B3 (Prolt+hardpt). We shall show that ha- 
los A and B have very different initial and final orbital properties 



despite the fact that their shapes in the presence of baryonic com- 
ponents are very similar (DOS). 



4.1 Distributions of orbital frequencies 

The frequency distribution of randomly selected orbits in a triax- 
ial halo can be used to characterise the orbital structure of phase 
space. It is useful to begin by discussing our expectations for how 
orbital frequencies change in response to growth of a central bary- 
onic component. The potential is significantly deeper in phase b 
compared to phase a, consequently the most tightly boimd orbits 
in the initial potential increase their frequencies. In contrast orbits 
which largely lie outside the central mass concentration do not ex- 
perience much deformation or much change in their frequencies. 
The higher the initial frequency the greater will be the frequency in- 
crease. Hence we expect a faster-than-linear increase in frequency 
in phase b relative to the frequency in phase a. 

When the baryonic component is evaporated, the halo expands 
once more and the halos regain their triaxiality in models SAl, 
Pi A3 and P/B2 but are irreversibly deformed in runs PiA4 and 
P;B3. One way to investigate the cause of the difference in these 
behaviours is to look for correlations between largest fundamen- 
tal frequencies of each orbit in each of the three phases. When the 
growth of the baryonic component causes an adiabatic change in 
orbits, one expects that their frequencies ui, will change in a reg- 
ular (i.e. monotonic) way so that the particles which are deepest 
in the potential experience the greatest frequency increase. In Fig- 
ures |4] and |5] we plot correlations between the frequencies uja, ij^t 
and LOc- 

In Figures|4]we show results for the three models whose bary- 
onic scale length is greater than 1 kpc: the left panels show that ut 
increases faster- than-linearly with cja, as expected with fairly small 
scatter. The right hand panels show that ujc is quite tightly corre- 
lated with LJa in all three models with the tightest correlation for 
simulation SAl (the dashed line shows the 1:1 correlation between 
the two frequencies). The deviation from the dashed line and the 
scatter is only slightly larger in simulations P!A3 and P/B2. The 
strong correlation between ojc and Ua in these models supports the 
argument by DOS that the growth of the baryonic component re- 
sulted in regular rather than chaotic evolution. In all three models 
only a small fraction of points deviate from the dashed line for the 
highest frequencies. 

In contrast, the two models with a hard central point mass, 
P;B3 and P; A4, (Fig. [Sj show a non-monotonic change in iot, in 
response to the growth of the central point mass as well as a higher 
degree of scattering in frequency space. In particular we note that 
orbits with small values of (i.e. those which are most weakly 
bound and have large apocenters) have the largest values of ujb, 
which is in strildng contrast to the situation in Figure|4] We also see 
that ujt sometimes decreased instead of increasing - again evidence 
for a scattering in frequency space rather than an adiabatic change. 
There is also greater scatter in the right-hand panels pointing to 
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Figure 4. For the three models with extended baryonic components the left panels show uj), versus uJa and right panels show uJc versus uia (frequencies in 
Gyr~^). Dashed lines in each panel show the 1:1 correlation between each pair of frequencies. From top to bottom the models contain a baryonic disk (SAl), 
a spherical bulge (P;A3) and a spherical elliptical (Py B2). 



a less complete recovery in the frequencies after the baryonic 
component is evaporated. Thus we see that when the central point 
mass is hard and compact there is significant orbit scattering. 

It is clear from a comparison of Figures |4]and [5] that a bary- 
onic component with a scale length of i?;, ~ 1 kpc or larger gener- 
ally causes a regular adiabatic change in the potential while a hard 
point mass (_R(, ~ 0.1 kpc) can produce significant chaotic scatter- 
ing. 

In both Figures |4]and |5]we see in the right hand panels that 
Uc < OJa especially at large values of coa (i-e. all points in the fig- 
ures lie systematically below the line, indicating a decrease in the 
frequency ujc)- Thus particles must have gained some energy im- 



plying that there has been a slight expansion in the DM distribution 
following the evaporation of the baryonic component. 

What fraction of orbits experience a large fractional change in 
frequency of an orbit from phase a to phase b, and from phase a to 
phase cl To investigate this we define: 

l^UJab = \{ljJa — i^b)/i^a\ (6) 
AcJac — \{ijJa — ljJc)/ljJa\- (V) 

The first quantity is a measure of the change in frequency distri- 
bution of orbits induced by the presence of the baryonic compo- 
nent, while the latter quantity measures the irreversibility of the 
evolution following the "evaporation of the baryonic component". 
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Figure 5. Same as Fig.|4] but for the two models with hard point masses of 0.1 kpc softening. Top panels show the effect in the triaxial halo (dominated by 
box orbits) and bottom panel shows the effect on the prolate halo (dominated by L-tube orbits). 




Figure 6. Kernel density distributions of Awaft {left panel) and of Awac {right panel) for particles in the halos SAl, P/B2, P;A3, PiA4, and PiB3 and as 
indicated by line legends. Each curve normalised to unit integral. The inset shows the full histogram for PiA4 plotted on a different scale. 



In Figure |6] we plot kernel density histograms of the distribution 
of the frequency change Aujab (left panel) and Acjac (right panel) 
for orbits in all five models as indicated by the line-legends. Each 
curve is normalized so that the area under it is unity. 



The distribution of Aujab is much wider for models PiA3, 
P; A4 and P;B3 than for the other two models. In these three models 
the scale-length of the baryonic component is ^ 1 kpc and results 
in a broad distribution of Auab, indicating that orbits over a wide 
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range of frequencies experience significant frequency change. For 
model P;A4 (dotted line) the histogram of values of AiJab appears 
almost flat on the scale of this figure because it is spread out over 
a much larger range of abscissa values indicating that many more 
particles are significantly scattered in phase b. The full distribution 
for PiA4 (see inset panel) is similar in form to P; A3 and PiB3. 

The right panels show that only a small number of orbits in 
models SAl, P;A3 and P/B2 experience an irreversible frequency 
change Auac > 20%, with the majority of particles experiencing 
less than 10%. In contrast in models P;A4 and P;B3, the distribu- 
tion of AiOac is much broader: a significant fraction of particles 
have experienced a large (20-50%) permanent change in their fre- 
quencies, reflecting the fact that the models with a hard-compact 
point mass are the only ones which experience irreversible chaotic 
scattering. 

Are there specific orbital characteristics that contribute to a 
large permanent frequency change, Aujac, between the two triaxial 
phases? We address this by determining how this quantity relates to 
other orbital properties. In Figure|7]we plot Auuac versus uja (left 
panels) and versus ujb (right panels) for four of our models. In the 
top two panels (SAl, and P/B2 - models with an extended bary- 
onic component) there is no evidence of a dependence of frequency 
change on and only a slight increase in Aujac at the highest val- 
ues of oja (results for P;A3 are not shown but are very similar to 
those for SAl). 

On the other hand, the lower two panels (P;A4 and PiB3 - 
models with a compact hard baryonic component) show that there 
is a strong correlation between Acjac and orbital frequency cua in- 
dicating that the orbits with the highest frequencies (cja) experi- 
ence the largest frequency change, Acjac- This is evidence that 
scattering by the hard central point mass is greatest for particles 
that are most tightly bound and therefore closest to the central po- 
tential, confirming prev ious expectations (Gerhard & Binnev 1985; 
iMerritt & Valluriiri996l) . The absence of an appreciable correlation 
with ujb is the consequence of scattering of orbits in frequency. 

In Figure[8]we plot Acjac versus rpcri (left panels) and versus 
litotl (the total specific angular momentum of an orbit averaged 
over its entire orbit in phase a) (right panels) for orbits in the two 
models with compact central point mass (P;A4 and PiB3). (We do 
not show plots for models SAl, P;A3 and P/B2, because they show 
no correlation between Acjac and either |jtot| or rperi.) The left 
panels of Figure [8] shows that orbits which pass closest to the cen- 
tral point mass experience the most significant scattering. The ab- 
sence of a correlation with | jtot j however indicates that scattering is 
independent of the angular momentum of the orbit. In the next sec- 
tion we show that halo A (the initial triaxial halo for model P;A4) 
is dominated by box orbits while halo B is initially prolate, and 
is dominated by L-tubes which circulate about the long axis. Con- 
trary to the prevailing view that centrophilic box orbits are more 
strongly scattered by a central point mass than centrophobic tube 
orbits, these figures provide striking evidence that chaotic scatter- 
ing is equally strong for the centrophobic L-tubes that dominate 



model PiB3 as it is for the centrophilic box orbits that dominate 
model P;A4. We return to a fuller discussion of the cause of this 
scattering in §|5] 



4.2 Changes in Orbital Classification 

As we discussed in § 13.21 relationships between the fundamental 
frequencies of a regular orbit can be used to classify it as a box 
orbit, a L-tube or a S-tube orbit. Quantifying the orbital composi- 
tion of the two different halos A and B and how their compositions 
change in response to the growth of a baryonic component yields 
further insight into the factors that lead to halo shape change. Or- 
bits were first classified as regular or chaotic based on their drift 
parameter log(A/) as described in §[3T| Regular orbits were then 
classified into each of three orbital families using the classification 
scheme outlined in § 13.21 The results of this orbit classification for 
each model, in each of the three phases, are given in Table|2] 

The most striking difference between the initial triaxial mod- 
els is that halo A (phase a of models SAl, PiA3 and PiA4) is dom- 
inated by box orbits (84-86%) while halo B (phase a of models 
P/B2 and P;B3) is dominated by L-tubes (78%). (The small differ- 
ences between models SAl, P!A3 and PiA4 in phase a is purely a 
consequence of the selection of different subsets of orbits from halo 
A.) None of the initial models has a significant fraction of S-tubes 
or of chaotic orbits. 

The very different orbit compositions of halos A and B in 
phase a results in rather different evolutions of their orbital pop- 
ulations in response to the growth of a central baryonic component. 
Although the growth of the disk results in a significant decrease 
in the box orbit fraction (from 86% to 43%) with boxes being con- 
verted to either S-tubes or becoming chaotic in phase b, model SAl 
is highly reversible suggesting an adiabatic change in the poten- 
tial. In model P;A3 and P; A4, the more compact spherical baryonic 
components decrease the fraction of box orbits even more dramat- 
ically (from 84% down to 16-17%), pointing to the vulnerability 
of box orbits to perturbation by a central component. Despite the 
similar changes in the orbital populations of the two models, P; A3 
is much more reversible than model P; A4, indicating that both the 
shape of the central potential and its compactness play a role in 
converting box orbits to other families and that the change in or- 
bit type is not evidence for chaotic scattering. It is striking that the 
more compact point mass in model P;A4 results in significantly 
more chaotic orbits (21%) compared to 8% in P;A3. 

While halo A is initial dominated by box orbits, halo B is ini- 
tially dominated by L-tubes, which dominate the orbit population 
in halo B in all three phases. The growth of the baryonic component 
in phase b causes the box orbit fraction to decrease (especially in 
model P!B3) while the fraction of chaotic orbits increases slightly. 
The more extended point mass in P/B2 causes a larger fraction of 
L-tubes to transform to orbits of another type than does the harder 
point mass in P!B3, despite the fact that there is much greater scat- 
tering in the latter model. A comparison between the model P;A4 
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Figure 7. For four models small dots show AoJac versus uJa (left) and versus lui, (right) for all particles analysed. Large solid dots with error bars show mean 
and standard deviation of particles in 15 bins in frequency. Top two panels are for models with an extended baryonic component (SAl, VfB2), lower two 
panels show models (P; A4, PiB3) with a compact baryonic component. 
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Table 2. Orbit composition of the models. The numbers represent the fraction of orbits in each family. 



Type 




Run SAl 






Run PiA3 






Run Pf 


B2 




Run P;A4 






RunPiB3 




Phase 


a 


h 


c 


a 


h 


c 


a 


b 


c 


a 


b 


c 


a 


b 


c 


Boxes 


0.86 


0.43 


0.83 


0.84 


0.16 


0.76 


0.15 


0.09 


0.29 


0.84 


0.17 


0.80 


0.15 


0.03 


0.21 


L-tubes 


0.11 


0.09 


0.12 


0.12 


0.43 


0.15 


0.78 


0.75 


0.54 


0.12 


0.35 


0.11 


0.78 


0.78 


0.59 


S-Tubes 


0.02 


0.27 


0.03 


0.02 


0.33 


0.06 


0.07 


0.09 


0.16 


0.02 


0.26 


0.04 


0.07 


0.11 


0.14 


Chaotic 


0.01 


0.21 


0.02 


0.02 


0.08 


0.03 


0.00 


0.07 


0.01 


0.02 


0.21 


0.05 


0.00 


0.08 


0.06 



and PiBS show that their orbit populations in the presence of a bary- 
onic component differ significantly due to the different original or- 
bit populations, while their degree of irreversibility is identical (e.g. 
Fig.|6} since the point mass in the two models is identical. 

A significant fraction (21%) of the orbits in phase b of model 
SAl and P;A4 are classified as chaotic (orbits with drift rate 



log (A/) > —1.0), in comparison with 9%, 7% and 8% in mod- 
els P;A3, P/B2 and PiB3 respectively. While the presence of such 
a large fraction of chaotic orbits in phase b of model P;A4 may 
be anticipated from previous work, the high fraction of chaotic or- 
bits in SAl (Triax+Disk) is puzzling. To address concerns about 
classification error that could arise from errors in the accuracy of 
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Figure 9. Distributions of rpcri for different orbit types. Distributions of each of the four different orbital types as indicated by the line-legends. Distribution 
in phase a is given by black curves and distribution in phase h is shown by red curves. The integral under each curve is proportional to the number of orbits of 
that orbital type. 



our frequency computation, we showed, in Figure [3] that chang- 
ing the frequency at which orbits were sampled by a factor of five 
did not result in any change in the overall distribution of log( A/), 
and hence should not affect our classification of orbits as regular or 
chaotic. Another puzzling fact is that, although model SAl mphase 
b has such a significant fraction of chaotic orbits, the orbit fractions 
essentially revert almost exactly to their original ratios once the 
disk is evaporated in phase c. Hence, the large fraction of chaotic 
orbits in phase b do not appear to cause much chaotic mixing. We 
will return to a more complete investigation of this issue in 14.41 

In Figure|9]we investigate how orbits of different types (boxes, 
L-tubes, S-tubes, chaotic) are distributed with rpcri, and how this 
distribution changes from phase a (black curves) to phase b (red 
curves). In phase a the initially triaxial halo A models (black 
curves) are dominated by box orbits. The fraction of box orbits 



is significantly decreased in phase b. In particular box orbits with 
large rpcri are transformed equally to short axis tubes and chaotic 
orbits, while some box orbits at small rpcri are converted to long- 
axis tubes. In contrast halo B models are dominated by L-tubes in 
both phases. Rather striking is how little the fraction of L-tubes 
in the halo B models changes, despite the fact that the halos are 
significantly more oblate axisymmetric in phase b than in phase 
a. We saw in Figure [6] that a significant fraction of orbits expe- 
rience strong scattering that manifests as a change in their orbital 
frequencies, and in Figure[8]we noted that the orbits with the small- 
est pericenter radii experience the largest change in frequency. In 
both models PiA4 and P;B3 the compact central point mass sig- 
nificantly reduced the box orbits. However model P;B3 only has 
a small fraction (15%) of box orbits and it seems unlikely that 
the chaotic scattering of this small fraction of orbits off the cen- 



© 2009 RAS, MNRAS Q0Q.mt24l 



The Orbital Evolution Induced by Baryonic Condensation in Triaxial Halos 1 5 



tral point mass is entirely responsible for driving the evolution of 
halo shape. Also P/B2 which has a much more extended baryonic 
component shows a change in orbit population which closely par- 
allels PjBS and we saw that P/B2 is quite reversible and shows 
little evidence for chaotic scattering. It is clear (from Fig.[9l( that in 
the prolate models (halo B) the majority of the orbits are L-tubes 
with large pericenter radii ((rperi) ~ 3 kpc) and these remain L- 
tubes in phase b. How then do these prolate models evolve to more 
spherical models while retaining their dominant orbit populations? 
To address this question we will now investigate the distribution of 
orbital shapes in each model at each phase of the evolution. 



4.3 Changes in orbital sliape 

A parameter, to quantify the shape of an orbit was defined in 
Equation 3 of § 13.31 Recall that this quantity is positive when the 
orbit is elongated along the major axis of the triaxial figure, is neg- 
ative when elongated along the intermediate axis, and almost zero 
when the orbit is "round" (ljx ^ Uy ^ ujz) or roughly axisym- 
metric about the minor axis {ux ~ ijJy)- In Figure [TO] we show 
the shape distributions for the orbits in four of our five models. For 
each model we show kernel density histograms for models in phase 
a (solid curves), phase b (dot-dashed curve), and phase c (dashed 
curves). In each plot the curves are normalized such that the inte- 
gral under each curve is unity. We define orbits to be elongated if 
Xa ^ 0.25, and to be "round" if |xs| ^ 0.1. 

Before the growth of the baryonic component (phase a: solid 
curves) the halo A models (left panels SAl, P;A4) have a distri- 
bution of orbital shapes that has a large peak at Xs ~ 0.35, arising 
from elongated orbits and a very small peak at Xs ~ due to round 
orbits (model P;A3 is not shown but is similar to P!A4.) In halo B 
models, (right panels P/B2, PiB3) on the other hand, the distribu- 
tion of shapes is double peaked with about one third of all orbits 
contributing to the peak at Xs ~ 0. This implies one third of its 
orbits in the initially prolate halo B are "round". In both halo A and 
B however, the larger of the two peaks has a value of Xs ~ 0.35 
corresponding to quite elongated orbits. Despite the quite different 
underlying orbital distributions (halo A models dominated by box 
orbits while the halo B models are dominated by L-tubes). This 
illustrates that despite having different orbital compositions, a sig- 
nificant fraction of their orbits are similarly elongated. 

The dot-dashed curves in all the panels show the distribution 
of orbital shapes in phase b. In all four models there is a dra- 
matic increase in the peak at jxs | ~ 0, pointing to a large increase 
in the fraction of round (or S-tubes) at the expense of the elon- 
gated (L-tube or box) orbits. In the halo B models the elongated 
orbits are significantly diminished indicating that the elongated L- 
tubes in phase a are easily deformed to "round" orbits in phase 
b (most likely squat inner-L-tubes). However, in model SAI there 
is a large fraction of orbits with intermediate values of elongation 

0.1 X. 0.4. 

In phase c (dashed curves) all models show the dominant peak 



shifting back to quite high elongation values of Xs ~ 0.3 (although 
this is slightly lower than Xs ~ 0.35 in phase a). The downward 
shift in the peak is most evident in model PiB3 (Prolt+hardpt), 
which as we saw before, exhibits the greatest irreversibility in 
shape. The scattering of a large fraction of the orbits by the hard 
central potential in model P;B3 seen in Figures|7]and[8]is the major 
factor limiting reversibility of the potential. The smallest shift is for 
model SAl (Triax+Disk), which exhibited the greatest reversibility. 

We can also investigate how the shapes of orbits vary with 
pericentric radius. We expect that orbits closer to the central poten- 
tial should become rounder (xs — >■ 0) than orbits further out. We 
see that this expectation is borne out in Figure [TT] where we plot 
orbital shape parameter Xs versus rperi in both phase a (left hand 
plots) and phase b (right hand plots). In each plot the dots show 
values for individual orbits. The solid curves show the mean of the 
distribution of points in each of 15 bins in rperi- Curves are only 
plotted if there are more than 30 particles in a particular orbital 
family (P;A4 is not shown since it is similar to PiA3). For models 
SAI (Triax+Disk) and PiA3 (Triax+Bulg) the figure confirms that 
elongated orbits in the initial halo A were box orbits (black dots 
and curves) and L-tubes (red dots and curves). The S-tubes (blue 
dots and curves) are primarily responsible for the "round" popula- 
tion at Xs ~ 0. In phase b (right-hand panels) of both SAl and 
P;A3 there is a clear tendency for the elongated orbits (boxes, L- 
tubes and chaotic) to become rounder at small pericenter distances, 
but they continue to be somewhat elongated at intermediate to large 
radii. Chaotic orbits in phase b of model SAI appear to span the full 
range of pericentric radii and are not confined to small radii. (Note 
that the density of dots of a given colour is indicative of the number 
of orbits of a given type but the relative fractions are better judged 
from Fig.|9]and Table 2.) 

For phase a in the models P/B2 (Prolt+Ellip) and P;B3 
{Prolt+hardpt) (left panels of each plot), boxes and L-tubes are 
elongated (Xs 5S 0.25), except at large rperi S5 8 kpc where they 
become rounder. We see a trend for the average orbital shape (as 
indicated by the curves) in phase b to become round at small peri- 
center radii. 

Note that in all the plots the curves only show the average 
shape of orbits of a given type at any radius. The points show that 
in the case of the L-tubes in particular, the red dots tend to be dis- 
tributed in two "clouds": one with large elongations Xs > 0.3 and 
one with small elongation Xs ~ 0.1. 

Thus in all four models it is clear that orbits that are elongated 
along the major axis of the triaxial potential in phase a become 
preferentially rounder at small pericenter radii in phase b. It is this 
change in orbital shape that plays the most significant role in caus- 
ing the overall change in the shape of the density in the baryonic 
phasfl 



Due to our chosen definition of shape pai'ameter, S-tubes generally have 
Xs regardless of radius, because uix ~ ujy. 
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Figure 10. Kernel density histograms of the distribution of orbital shape parameter Xs for each of the four models: SAl (top left), P;A4 bottom left , P/B2 
(top right ) and P;B3 (bottom right) (P; A3 is not shown since it is very similar to P;A4.) Distributions of Xb in phase a are shown by solid curves, in phase b 
by dot-dashed curves, and in phase c by dashed curves. In all models, a large fraction of orbits in phase b are "round" (Xs — 0..) 



4.4 Frequency maps and chaotic orbits 

We saw in Table |2] that phase b of model SAl (Triax+Disk) and 
of model PiA4 (Triax+hardpt) have a significant fraction (21%) 
of chaotic orbits (i.e. orbits with log(A/) > —1). While P;A4 
shows significant lack of reversibility, which we can attribute to the 
presence of this high fraction of chaotic orbits, model SAl does not 
show evidence for irreversibility. 

Figure[72lshows kernel density histograms of the chaotic drift 
parameter log(A/) for orbits in each of the three phases in model 
SAl. It is obvious that in phases a and c there is only a small frac- 
tion of chaotic orbits (i.e. orbits with log(A/) > —1), whereas a 
much more significant fraction of orbits lie to the right of this value 
in phase b. Even the peak of the distribution in phase b is quite 
significantly shifted to higher drift values. 

In this section we investigate the surprising evidence that the 



chaotic orbits in phase b of model SAl do not appear to mix. One 
possible reason for the lack of diffusion of the chaotic orbits is that 
the timescale for evolution is not long enough. Indeed DOS report 
that evolving run SAl with the disk at full mass for an additional 
5 Gyr after the growth of the disk is complete, leads to a larger 
irreversible evolution (see their Figure 3a). Nonetheless, even in 
that case the irreversible evolution was only marginally larger than 
when the disk was evaporated right after it grew to full mass. More- 
over the growth time was 5 Gyr which means that the halo was 
exposed to a massive disk for a cosmologically long time. 

A second possible reason for the lack of chaotic diffusion is 
that most of the chaotic orbits in this phase of the simulation are 
"sticky". The properties of "sticky chaotic orbits" were discussed 
in §[3] In a series of experiments designed to measure the rate of 
chaotic mixing, Merritt & Valluri (.1996.) showed that while ensem- 
bles of strongly chaotic orbits diffused and filled an equipotential 
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Figure 11. For models SAl, PjAS, P/B2 and PjBS, the orbital shape parameter Xs for each orbit is plotted against its pericentric radius rpcii as a small dot. 
The orbits of each of the four major orbital families are colour coded as in the figure legends. Left hand panels are for phase a and right hand panels are for 
phase b. The solid curves show the mean value of Xs for all particles of that particular family, in 15 bins in Tpori- Curves are not plotted if there are fewer than 
30 orbits in a given orbital family. We used a kernel regression algorithm to smooth the cui'ves. 
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Figure 12. Histograms of frequency drift parameter log(A/) for tlie three 
phases of model SAl as indicated by the line-legends. 

surface on timescales between 30-100 dynamical times, similar en- 
sembles of "sticky" or resonantly trapped orbits diffused much less 
quickly and only filled a small fraction of the allowed surface after 
very long times. 

Laskar (1990) showed that frequency maps are a powerful way 
to identify resonances in dynamical systems. Frequency maps are 
obtained by plotting ratios of the 3 fundamental frequencies for 
each individual orbit. If a large and representative orbit population 
is selected, they can provide a map of the phase space structure 
of the potential including all the resonances. Resonances appear 
as straight lines on the frequency map since their fundamental fre- 
quencies satisfy a condition like luj^ + mujy + noj^ = 0. This 
method of mapping the phase space has the advantage that since it 
only depends on the ratios of the frequencies and not on the fre- 
quencies themselves, it can be used to map phase space for large 
ensembles of particles without requiring them to be iso-energetic. 
This is a significant advantage over mapping schemes like Poincare 
surfaces-of-section, when applied to an TV-body simulation where 
particles, by design, are initialised to be smoothly distributed in en- 
ergy. Thus one can use the method to identify global resonances 
spanning a large range of orbital energies in TV-body simulations. 

In Figure [T3] we present frequency maps for four of the five 
models in phase a (left panels) and phase b (right panels) (P/B2 
is not shown since the frequency maps for this model are indistin- 
guishable from those for P;B3). For each orbit the ratios of the fun- 
damental frequencies ujy/uj^ and uj^/^j^z are plotted against each 
other. Particles are colour coded by their energy in phase a. The 
energy range in phase a was divided into three broad energy bins, 
with equal numbers of particles per bin. The most tightly bound 
particles are coloured blue, the least bound particles are coloured 
red and the intermediate energy range is coloured green. 

Resonance lines are seen in the clustering of particles in all the 
maps. The most striking of the frequency maps is that for phase b of 



model SAI (Triax+Disk). This map has significantly more promi- 
nent resonance lines, around which many points cluster, than any 
of the other maps. Three strong resonances and several weak reso- 
nances are clearly seen as prominent straight lines. The horizontal 
line al ujy/ijjz ~ 1 corresponds to the family of orbits associated 
with the 1:1 closed (planar) orbit that circulates around the a;-axis, 
namely the family of "thin shell" L-tubes. The diagonal line run- 
ning from the bottom left corner to the top right comer with a slope 
of unity (ujy/uuz = i^i/i^z) corresponds to the family of orbits 
that circulates about the z-axis: the family engendered by the "thin 
shell" S-tubes. Since this latter family shares the symmetry axis of 
the disk, it is significantly strengthened in phase b by the growth 
of the disk. In addition to having many more orbits associated with 
it, this resonance extends over a much wider range in energy as 
evidenced by the color segregation along the resonance line (blue 
points to the bottom left and red points at the top right). This seg- 
regation is the result of an increase in the gradient of the potential 
along the z-axis due to the growth of the disk, which results in 
an increase in lj^- The more tightly bound a particle, the greater 
the increase in ljz, and the greater the decrease in both its ordi- 
nate and abscissa. The most bound particles (blue points) therefore 
move away from their original positions towards the bottom left 
hand corner of the plot. The least bound particles (red points) are 
furthest from the center of the potential and these points experience 
the least displacement - although these points also shift slightly to- 
ward the resonance lines. 

A third prominent resonance is the vertical line at lo^l^z = 
0.5 that corresponds to orbits associated with the family of banana 
(1:2 resonant) orbits. This banana (boxlet) resonance is also en- 
hanced by the growth of the disk since this family of orbits, while 
not axisymmetric, is characterised by large excursions along the 
X-axis and smaller excursions in the z direction. Several shorter 
resonance lines are seen but are too sparsely populated in this plot 
to properly identify. 

The frequency map for model PiA3 in phase b shows that 
a spherical baryonic component produces a rather different phase 
space structure than that produced by the disk. In particular it is 
striking that the most tightly bound (blue) points are now clustered 
at the intersection of the horizontal and diagonal resonances namely 
around the closed period orbits 1:1:1. This may be understood as 
the consequence of the growth of the spherical baryonic point mass 
around which all orbits are rosettes and since no direction is pre- 
ferred all orbits are "round". The 1:2 banana resonance is also less 
prominent in this model (largely because the deep central potential 
destabilises this boxlet family). 

The frequency map for model P; A4 phase b shows the greatest 
degree of scattering, as evidenced by the thickest resonance lines. 
We attribute this to the large number of chaotic orbits in this model. 
Apart from the broad clustering of points around the diagonal (S- 
tube) and horizontal (L-tube) resonances there are no strong reso- 
nance lines seen in this map. Unlike the map for P; A3 which shows 
a clustering of tightly bound (blue) points at the 1:1:1 periodic or- 
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Figure 13. Frequency maps of particles in phase a and phase h for four models. For each particle the ratio of the fundamental frequencies uiy/uiz is plotted 
versus lo^I^z is plotted by a single dot. The dots are colour coded by the energy of the particle in phase a. The most tightly bound particles are coloured blue, 
and the least bound particles are coloured red. Model SAl has 6000 particles, model P; A4 and PjAS have 5000 particles, while PjES has 1000 particles each. 



bit resonances, the blue points are widely scattered in the frequency 
map of P;A4. 

The frequency maps for model PiB3 shows that most of the 
orbits in this model are associated with the (1:1) L-tube family 
(horizontal line). A smaller number of orbits is associated with the 
1:1 S-tube resonance (diagonal line). We saw previously that the 
growth of the baryonic components in this prolate halo caused lit- 
tle change in the orbit families. This is confirmed by the fact that the 
frequency maps in both phases are remarkably similar except for an 
increase in the clustering of points at the intersection of the hori- 
zontal (L-tube) and the diagonal (S-tube) resonance, which occurs 
for the same reason as in PiA3. Since halo B is initially a highly 
prolate model, it has (as we saw previously) only a small fraction 
of box (and boxlet) orbits and in particular no banana orbits. 

It is quite striking that in phase a the frequency maps show 
significantly less segregation by energy, and only a few resonances. 
This is because the initial triaxial models were generated out of 
mergers of spherical NFW halos which were initially constructed 
so that orbits were smoothly distributed in phase space. The in- 
crease in the number of resonances following the growth of a bary- 
onic component is one of the anticipated consequences of resonant 
trapping that occurs during the adiabatic change in a potential (e.g. 
[rremaine & Yu 2000; Binney & Tremaine 2008). 

To test the conjecture that the majority of chaotic orbits in 



model SAl (phase b) are resonantly trapped, we compute the num- 
ber of chaotic orbits that lie close to a major resonance line. We de- 
fine "closeness" to the resonance by identifying those orbits whose 
frequency ratios lie ±a of the resonant frequency ratio. For exam- 
ple we consider an orbit to be close to the (1:1) L-tube resonance 
(horizontal line in map), if \uiy/u}z — 1| ^ a. We find that the frac- 
tion of chaotic orbits in phase b, that lie close to one of the three 
major resonances identified above, is 5 1% when a — 0.01 and 62% 
when Q — 0.03. Weaker resonances lines (which are hard to recog- 
nise due to the sparseness of the data points) may also trap some of 
the chaotic orbits. This supports our conjecture that the main reason 
that model SAl does not evolve in phase b, despite the presence of 
a significant fraction of chaotic orbits, is that the majority of the 
chaotic orbits are trapped around resonances and therefore behave 
like regular orbits for very long times. 

In Figure[T4lwe plot four examples of chaotic orbits in phase b 
of SAl, which illustrate how resonantly trapped or "sticky" chaotic 
orbits look. Each panel of four sub-plots shows a single orbit plot- 
ted in two Cartesian projections (side-by-side). The top pair of sub- 
plots show the orbit over the first 10 Gyr long time segment, while 
the bottom pair shows the same orbit over a second 10 Gyr time 
segment. The two time segments were separated by 10 Gyr. For il- 
lustration we selected orbits with a range of drift parameters. The 
orbit in the top-left panel is an example of an orbit that conforms to 
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Figure 14. Several chaotic orbits in phase h of model SAl. Orbits are plotted in two Cartesian projections over two different time segments of 10 Gyr (top two 
plots in each panel show the first time segment, while bottom panels show the second time segment). See text for details. 



our notion of a chaotic orbit that explores more phase space as time 
progresses, and has a large drift parameter of log(A/) = —0.48. 
The top-right panel shows a S-tube orbit that suddenly migrates to 
a box orbit (this orbit has a log(A/) = —0.54) and was proba- 
bly in the separatrix region between the S-tube and box families. 
The bottom-left panel shows an orbit that is originally a S-tube that 
becomes trapped around a resonant boxlet ("fish") family (with 
log(A/) = —0.66), while the bottom right-hand panel shows a 
weakly chaotic box orbit (with log(A/) = -0.94). Of the 21% 
of orbits in phase b that are chaotic (log(A/) ^ —1), only ~ 5% 
have log(A/) ^ —0.5. This fraction is small enough that one does 
not expect it to result in significant chaotic mixing. 



5 SUMMARY AND DISCUSSION 

Since it was first proposed, the idea that a central black hole would 
scatter centrophilic box orbits in triaxial galaxies resulting in more 
axisymmetric potentials (Gerhard & Binnev 1985) has frequently 
been used to explain the shape change in a va riety of systems from 
the destruction of bars by central black holes jNorman et al.lll99^ 
to t he formation of more obl ate galaxy clusters in simulations with 
gas jKazantzidis et alj2004l) . 

Experiments on chaotic mixing indicated that the 
timescales for such mi xing is about 30 - 100 dynamical times 
jMerritt & Vallurill996h . which is much longer than the timescales 
for evolution of dark matter halos in simulations with gas 
jKazantzidis et alj |2004|) . In addition recent detailed studies of 
A'^-body simulations with controlled experiments have shown that 
the role of chaotic mixing may be less dramatic than conjectured 
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by these previous studies. A study of relaxation of collisionless 
systems following the merger of two spherical galaxies showed 
that despite the fact that a large fraction of the orbits in a system 
undergoing violent relaxation are chaotic, the timescales for 
chaotic diffusion and mixing are too long for this process to play 
a significant role iValluri et alj|2007h . In fact, even after violent 
relaxation, orbits retain strong memories of their initial energies 
and angular momenta. 

DOS argued that since chaos introduces an irreversible mix- 
ing, numerical experiments in which evolution is driven by chaotic 
orbits should not be reversible. These authors studied the macro- 
scopic shapes of triaxial dark matter halos in response to the growth 
of a baryonic component. Unless the baryons were too centrally 
concentrated, or transported angular momentum to the halo, the 
evolution they saw was reversible, from which they concluded that 
much of the shape change arises from deformations in the shape 
of individual orbits rather than significant chaotic scattering. In 
this paper we investigated this issue in significantly greater de- 
tail by applying the Numerical Analysis of Fundamental Frequen- 
cies (NAFF) technique that allows us to quantify the degree to 
which chaotic diffusion drives evolution and to identify the pri- 
mary physical processes that cause halo shape change. The fre- 
quency based method is able to distinguish between regular and 
chaotic orbits, making it more useful than Lyapunov exponents 
which ar e known to be sensitive to dis cretization effects in A'^-body 
systems l lHemsendorf & Merrittl20o3) . We use the method to quan- 
tify the drift in frequencies of large representative samples of orbits, 
thereby quantifying the degree of chaos in the systems we study. It 
also allowed us to map the phase space structure of the initial and 
final halos and to quantify the relationship between the change in 
the shapes of individual orbits and the shape of the halo as a whole. 

Applying various analysis methods to orbits in five systems 
we demonstrated that the conclusion reached by DOS that chaos is 
not an important driver of shape evolution when the baryonic com- 
ponent is extended is indeed valid. As did DOS, we also found that 
significant chaotic scattering does occur when the baryonic com- 
ponent is in the form of a hard central point mass (of scale length 
^ 0.1 kpc). It is interesting that regardless of the original orbital 
composition of the triaxial or prolate halo, and regardless of the 
shape of radial scale length of the baryonic component, halos be- 
come more oblate following the growth of a baryonic component. 
Thus two quite different processes (chaotic scattering and adiabatic 
deformation) result in similar final halo shapes even in halos with 
very different orbital compositions. 

We explored two different initial halos, one in which box or- 
bits were the dominant elongated population (halo A) and the other 
in which L-tubes dominated the initial halo (halo B). Despite the 
different orbit compositions both models exhibit similar overall 
evolution with regard to the shapes of orbits. In the halo A models, 
the box orbits were much more likely to change to either L-axis 
or S-tubes, whereas in the halo B models, the dominant family of 



L-tubes largely retained their orbital classification while deforming 
their shapes. 

Below we list the main results of this paper: 

(i) Correlations between the orbital frequencies in the three dif- 
ferent phases uja, oJt and Uc are a useful way to search for regular 
versus chaotic evolution of orbits. The orbital frequencies in the 
three phases are found to be strongly correlated with each other 
when the baryonic component is extended (Fig. |4), but show sig- 
nificant scattering when the baryonic component is a compact point 
mass (of scale length about 0.1 kpc) (Fig.O. In the more extended 
distributions, only a small fraction of the orbits experience signifi- 
cant change in their original orbital frequencies when the baryons 
are evaporated, while both the magnitude of scattering in orbital 
frequency as well as the fraction of orbits experiencing scatter- 
ing, increases as the baryonic component becomes more compact 
(Fig.[6j. 

(ii) In the three models with relatively extended baryonic com- 
ponents, the change in orbital frequency between phase c and phase 
a (AoJac) is not correlated with orbital frequency (Fig. |7), peri- 
center distance or orbital angular momentum. When the baryonic 
component is a hard point mass, however, the frequency change 
is greater for orbits that are deeper in the potential and therefore 
have both a higher initial orbital frequency (Fig. |7j and smaller 
pericentric radius (Fig.[8). Scattering in frequency affects both the 
centrophilic box orbits as well as centrophobic L-tubes. 

(iii) The growth of a baryonic component in halo A (either disk 
or softened point mass) causes box orbits with large pericentre radii 
to be converted to S-tubes, L-tubes or become chaotic (Fig|9] top 
panel). While this change is almost completely reversible in the 
case of the disk or a diffuse point mass, it is less so when the bary- 
onic component is a hard point mass. In halo B, which is dominated 
by L-tubes, the growth of the baryonic component causes almost no 
change in the orbital composition of the halo, indicating that the L- 
tubes are not destroyed but deformed (Fig|9] bottom panel). Even 
though P;B3 (Pwlt+hardpt) is a model with significant orbit scat- 
tering by the hard central point mass, the process appears to mainly 
convert elongated inner L-tube orbits to somewhat rounder outer 
L-tubes. In model P; A4 (Triax+hardpt) the box orbits are scattered 
onto S-tubes or chaotic orbits. The significant amount of scattering 
seen for even centrophobic L-tube orbits shows that the evolution 
is not due to direct scattering by a central point mass as sometimes 
assumed. Two alternative possibilities are more likely to account 
for the significant scattering in frequency. First the change in the 
symmetry and depth of the central potential is a perturbation to 
the potential that gives rise to an increase in the region of phase 
space occupied by resonances (Kandrup 1998 - private communi- 
cation). As the reson ances overlap th ere is an increase in the degree 
of chaotic behavior ( IChirikov|[T979l) . The second option is that the 
point mass attains equipartition with the backgrou nd mass distribu- 
tion, resulting in Brownian motion jMerrit t' 2005). The Brownian 
motion can cause the centre of the point mass to wander within a 
region of radius ~ 0.1 — 1 kpc which can result in a significant 
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change in the maximum gravitational force experienced by an or- 
bit from one pericentre encounter to the next. This change in the 
maximum gravitational force manifests as scattering of the orbit 
which is equally effective for both box orbits and long-axis tubes. 
Indeed a small wandering of the central massive point is seen in 
the A'^-body simulations of P; A4; this motion is not included in our 
orbit calculations since all particles are frozen in place when calcu- 
lating orbits. While it is beyond the scope of this paper to explore 
this issue further, we caution that the motion of the point mass in 
our A'^-body simulations is likely to over-estimate the magnitude 
of the Brownian motion, since this depends on the mass resolution 
of the background particles. This suggests that in a real galaxy the 
evolution of the shape is much more likely to be driven by smooth 
adiabatic deformation of orbits than chaotic scattering. 

(iv) In triaxial halos, the orbital shapes sharply peaked distribu- 
tion with the most elongated orbits (Xs > 0.25) are either boxes or 
L-tubes. In the prolate halos the second peak at Xs ~ contains 
a third of the orbits and is composed of squat outer L-tubes and 
some box orbits. The growth of a baryonic component of any kind 
causes orbits of all types to become "rounder", especially at small 
pericenter radii. This change in orbital shape distribution with ra- 
dius is the primary cause of the change of halo shapes in response 
to the growth of a baryonic component. This is consistent with the 
findings of DOS who also found that the orbits in the models be- 
came quite round. 

(v) The growth of a disk causes a large fraction of halo orbits 
to become resonantly trapped around major resonances. The three 
most important resonances are those associated with the 1:1 tube 
(thin shell) orbit that circulates about the short axis in the x ~ y- 
plane, the 1 : 1 tube (thin shell) orbit that circulates about the long 
axis in the y — z-plane and the 1:2 banana resonance in the x — z- 
plane. We saw from the frequency maps that the resonant trapping 
of the halo particles depends both on the form of the baryonic com- 
ponent grown in the halo as well as on the initial orbital population 
of the halo. 

Thus we conclude that the evolution of galaxy and halo shapes 
following the growth of a central component occurs primarily due 
to regular adiabatic deformation of orbital shapes in response to the 
changing central potential. Chaotic scattering of orbits may be im- 
portant particularly for orbits with small pericentre radii but only 
when the central point mass is extremely compact. Contrary to pre- 
vious expectations, chaotic scattering is only slightly more effec- 
tive for centrophilic box orbits than it is for centrophobic L-tubes. 
Boxes can be scattered onto both L- and S-tube orbits and a sig- 
nificant fraction become chaotic. When the compact central point 
mass scatters L-tubes as it does in a prolate halo, they are scattered 
onto other L-tube orbits rather than onto S-tube orbits. The strong 
chaotic scattering that we see on centrophobic L-tube orbits has not 
been previously anticipated. 

An important implication of our analysis is that while the 
shapes of halos (and by extension elliptical galaxies) become more 
oblate (especially at small radii), following the growth of a baryonic 



component, the majority of their orbits are not S-tubes as might 
be predicted from their shapes. Instead our analysis shows that or- 
bits prefer to maintain their orbital characteristics, and the majority 
of the orbits are those which would be generally found in triaxial 
galaxies. This is particularly important for studies of the internal 
dynamics of elliptical galaxies since the fact that their shapes ap- 
pear nearly axisymmetric need not imply that their orbital structure 
is as simple as the structure of oblate elliptical galaxies. Modifying 
the shapes to slightly triaxial could result in significant changes in 
their orbit populations and consequently could affect both the in- 
ferred dynamical structure as well as the estimates of the masses 
components such as the supermass ive black holes in these galaxies 
Jvan den Bosch & de Zeeuwll2009h . 

Finally, our finding that the growth of a stellar disk can re- 
sult in a large fraction of halo orbits becoming trapped in reso- 
nances could have important implications for observational studies 
of the Milky Way's stellar halo. The computation expense of the 
orbit calculations forced us to restrict the size of the frequency map 
for model SAl to 6000 particles. This is only a tiny fraction of the 
particles in the original simulation. Despite the smallness of the 
sample, the frequency maps (Fig. [13) shows a rich resonant struc- 
ture which implies that the particles (either stars or dark matter) in 
the stellar and dark matter halos of our Galaxy, particularly those 
close to the plane of the disk, are likely to be associated with reso- 
nances, rather than being smoothly distributed in phase space (this 
is in addition to structures arising due to tidal destruction of dwarf 
satellites). Although significantly greater resolution is required to 
resolve such resonances than is currently available, this could have 
significant implications for detection of structures in current and 
upco ming surveys of the Mi l ky Way such as SDSS-III (Segue) and 
Gaia IPerrvman et al.ll200ll ; IwiUdnson & et al.l 120051 and in on- 
going direct detection experiments which search for dark matter 
candidates. 
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